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Abstract 


In  this  report,  we  present  preliminary  results  on  the  design  of  the  baseline  controller  for  an  F- 
16  aircraft  automatic  landing  system  using  linear  matrix  inequalities  (LMI)-based  ap¬ 
proaches.  We  start  with  a  general  study  of  aircraft  control  and  dynamics  to  gain  knowledge  of 
the  structure  of  an  aircraft  dynamic  model  and  its  inner  loop  control  system.  We  then  identify 
a  linear  model  along  the  glide  path  for  the  inner  loop  control  system  in  the  simulator.  With 
this  linear  model,  the  control  objective  is  to  solve  a  stabilization  problem — stabilizing  the 
aircraft  along  the  glide  path  using  linear  state  feedback  controls.  Expressing  the  stability  cri¬ 
terion  and  the  constraints  in  LMIs,  we  cast  the  stabilization  problem  as  an  optimization 
problem.  Using  the  SDPSOL1  software  package  developed  by  Wu  and  Boyd,  we  solve  this 
optimization  problem  for  the  control  gain  and  stability  region,  which  completes  the  controller 
design  [Wu  96]. 


1  SDPSOL  is  a  parser/solver  for  semidefinite  programming  and  determinant  maximization  problems 
with  matrix  structure. 
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1  Introduction 


As  part  of  the  Evolutionary  Design  of  Complex  Software  (EDSC)  Incremental  Software 
Evolution  for  Real-Time  Systems  (INSERT)  Project  [INSERT  00],  the  Simplex™  architec¬ 
ture  has  been  implemented  in  a  legacy  F-16  manned  virtual  desktop  simulation  to  demon¬ 
strate  its  capability  of  supporting  the  safe  insertion  of  new  technology  and  the  safe  upgrading 
of  current  functionality  in  safety-critical  real-time  systems.  Using  the  Simplex  architecture, 
new  capabilities  can  be  rapidly  installed  in  the  aircraft  to  allow  the  fighter  to  perform  with 
enhanced  maneuverability.  By  tolerating  possible  faults  that  such  installations  may  bring  to 
the  system,  the  Simplex  architecture  guarantees  continuous  operation  of  the  aircraft  even  if 
the  newly  configured  aircraft  has  not  been  fully  tested  in  its  operating  environment. 

Fault  tolerance  in  the  Simplex  architecture  is  based  on  the  concept  of  analytic  redundancy. 

The  Simplex  architecture  consists  of  multiple  “replacement  units”  where  different  technolo¬ 
gies  can  be  implemented  and  where  at  least  one  of  them  will  be  highly  reliable  in  the  sense 
that  it  can  serve  as  the  fallback  option  when  others  fail  to  function  properly.  (For  more  infor¬ 
mation  on  the  Simplex  architecture,  refer  to  [Feiler  99],  [Poliak  98],  and  [Altman  99].)  Where 
control  of  the  aircraft  is  concerned,  there  must  be  a  fallback  controller  in  the  system  to  sup¬ 
port  the  upgraded  control.  This  controller  will  be  referred  as  the  baseline  controller.  Since  the 
baseline  controller  takes  over  control  upon  detection  of  a  fault  caused  by  the  upgraded  con¬ 
troller,  it  is  important  that  the  aircraft  be  in  a  state  from  which  the  baseline  controller  is  capa¬ 
ble  of  maintaining  the  aircraft’s  stability.  Therefore,  in  addition  to  the  design  of  the  control 
algorithm  for  the  baseline  controller,  an  operational  region  of  the  controller  in  the  state  space 
of  the  aircraft  needs  to  be  identified.  An  operational  region  of  a  controller  means  that,  from 
any  state  inside  the  region,  the  controller  achieves  the  control  objectives  without  violating 
any  constraints  imposed  upon  the  aircraft.  Furthermore,  such  a  region  can  serve  as  a  switch¬ 
ing  criterion  for  fault  detection. 

In  this  report,  we  address  the  issues  of  control  algorithm  design  and  operational  region  identi¬ 
fication  for  a  baseline  controller.  In  particular,  we  concentrate  on  the  landing  operation.  We 
present  the  development  and  preliminary  results  of  the  design  of  a  baseline  controller  for 
automatic  landing  of  an  F-16  using  modem  control  techniques  and  linear  matrix  inequality 
(LMI)-based  approaches.  The  control  objective  is  to  maintain  the  aircraft  flying  along  the 
glide  path,  which  is  equivalent  to  a  problem  of  stabilization  at  an  aircraft  steady  state.  A  lin¬ 
ear  model  is  identified  using  data  of  typical  landing  maneuvers  from  the  legacy  F-16  simula¬ 
tion.  LMIs  are  used  to  generate  a  feedback  controller  taking  into  account  actuator  constraints 
and  flight  specifications.  As  a  result  of  the  stabilization  problem,  a  stability  region  of  the 

™  Simplex  is  a  trademark  of  Carnegie  Mellon  University. 
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closed-loop  system  is  identified  to  serve  as  the  operational  region  of  the  baseline  controller,  a 
stabilizer. 

The  work  reported  here  was  carried  out  by  using  data  from  the  F-16  Block  50  model  simula¬ 
tion,  referred  as  the  real-time  simulator  (or  the  simulator)  in  this  report,  in  parallel  with  a 
study  of  general  aircraft  control  introduced  by  Stevens  and  Lewis  and  the  F-16  modeling  and 
simulation  tools  provided  by  Stevens  [Stevens  92].  A  software  package  implementing  soft¬ 
ware  illustrations  in  Aircraft  Control  and  Simulation  by  B.  Stevens  was  obtained  from  the 
author.  This  package  contains  a  nonlinear  model  of  an  F-16,  with  which  a  linearized  model 
can  be  derived  and  nonlinear  simulation  can  be  performed.  To  distinguish  from  the  real-time 
simulator  and  the  model  within,  we  call  this  simulation  a  book  simulator  and  the  nonlinear 
model  a  book  model. 

Understanding  the  book  model  of  the  F-16  and  inner  loop  control  design  as  described  by 
Stevens  provided  us  a  structure  for  the  inner  loop  F-16  control  system  in  the  simulator,  which 
is  not  available  to  us  [Stevens  92],  With  this  structure,  we  were  able  to  identify  the  control 
system  in  the  simulator  to  which  a  control  algorithm  can  be  designed  and  implemented.  Fur¬ 
thermore,  with  the  book  model,  stability  analysis,  and  LMI-based  control  design  has  been 
completed  at  an  advanced  stage. 

This  report  is  organized  as  the  follows.  Section  2  gives  an  introduction  to  aircraft  dynamics 
and  the  notation  used  throughout  the  report.  Section  3  presents  the  problem  to  be  addressed. 
System  identification  issues  are  discussed  in  Section  4.  Section  5  describes  the  design  of 
control  algorithms  and  derivation  of  stability  regions.  A  summary  of  the  work  and  main  con¬ 
clusions  are  presented  in  Section  6. 
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2  Dynamics  and  Control  of  an  Aircraft 


In  this  section,  we  present  a  review  of  aircraft  dynamics  and  control  described  by  Stevens 
[Stevens  92],  We  first  introduce  the  state  variables  and  control  inputs  in  a  typical  aircraft  and 
various  coordinate  frames  used  to  describe  aircraft  dynamics.  Then  we  establish  sets  of  equa¬ 
tions  of  motion  for  an  aircraft  in  different  coordinate  frames  and  flight  conditions.  Finally,  we 
describe  a  software  package  for  F-16  model  trimming  and  linearization. 


2.1  Aircraft  State  Variables  and  Control  Inputs 

As  a  rigid  body  moving  in  a  three-dimensional  space,  an  aircraft  has  a  total  of  6  degrees  of 
freedom  described  by  12  state  variables.  These  variables  are  divided  into  four  groups: 

1.  three  position  variables  [position  of  the  aircraft  center  of  gravity  (eg)] 

2.  three  linear  velocity  variables  (translational  velocity  of  the  aircraft  eg) 

3.  three  attitude  (orientation)  variables 

4.  three  angular  velocity  variables 

Aircraft  dynamics  are  normally  controlled  by  four  physical  inputs:  throttle,  aileron,  elevator, 
and  rudder.  The  throttle  controls  the  thrust  to  the  aircraft,  while  the  aileron,  elevator,  and  rud¬ 
der  deflections  generate  aerodynamic  forces.  Figure  1  shows  a  sketch  of  an  aircraft  and  its 
deflector  components  as  well  as  the  moments  that  they  generate. 


Directional  control 
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Figure  1:  Deflection  Actuators  and  Moments 
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2.2  Coordinate  Frames 

To  describe  the  dynamics  of  an  aircraft,  we  need  to  set  up  a  coordinate  frame.  Depending  on 
the  flight  conditions,  different  coordinate  frames  may  be  used  to  describe  an  aircraft’s  mo¬ 
tion.  The  following  summarizes  some  of  the  coordinate  frames.  Refer  to  Figure  3  for  a  pic¬ 
ture  of  the  relationship  between  the  coordinate  systems. 

Earth-centered  inertial  (ECI)  reference  frame:  centered  at  the  origin  of  the  Earth. 

•  x-axis  along  the  Earth’s  spin  axis,  pointing  to  the  North  Pole 

•  y-axis  perpendicular  to  x-z  plane 

•  z-axis  pointing  to  the  1 80°  longitude  point  on  the  equator  at  t  =  0 

North-east-down  (NED)  frame:  centered  on  the  Earth’s  surface  at  the  point  vertically  below 
the  aircraft  eg,  with  x-y  plane  tangent  to  the  Earth’s  surface.  It  moves  with  the  aircraft. 

•  x-axis  pointing  to  the  north 

•  y-axis  pointing  to  the  east 

•  z-axis  normal  to  the  Earth’s  surface  and  pointing  inward 

Aircraft-body  coordinate  (ABC)  frames:  refers  to  three  coordinate  frames,  which  are  all 
centered  at  the  aircraft  eg  with  the  conventional  right-handed  (forward,  starboard,  and  down) 
set  of  axes.  These  frames  are  illustrated  in  Figure  2. 

•  Body-fixed  axes  (ABC):  All  axes  are  fixed  on  the  body  of  the  aircraft. 

•  Stability  axes  (ABC-Stability):  Axes  are  obtained  from  the  body-fixed  axes  by  rotating 
about  the  y-axis  for  the  angle  of  a ,  which  is  called  the  angle  of  attack.  It  is  positive  if  the 
rotation  about  the  body-fixed  y-axis  is  negative. 

•  Wind  axes  (ABC- Wind):  Axes  are  obtained  from  the  stability  axes  by  rotating  about  the 
z-axis  for  the  angle  of  /?  ,  which  is  called  the  side-slip  angle.  It  is  positive  if  the  rotation 
about  the  stability  z-axis  is  positive. 

The  angles  a  and  /?  are  known  as  the  aerodynamic  angles;  these  are  needed  to  specify  the 
aerodynamic  forces  and  moments. 

Different  coordinate  frames  will  be  used  to  simplify  the  description  of  the  aircraft  dynamics 
simple  and  to  make  them  easier  to  understand.  For  example,  to  describe  a  flight  around  the 
Earth,  the  ECI  frame  is  usually  used,  while  flying  on  a  flat  earth  is  often  described  in  the 
NED  frame.  Algebraic  relations  among  these  frames,  namely  the  mapping  from  one  frame  to 
another,  are  given  in  Appendix  A.  Figure  3  shows  the  spatial  relationship  among  the  ECI 
frame,  NED  frame,  and  ABC  frame. 
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Figure  2:  Aircraft  Body  Coordinate  System 

y 


Figure  3:  Relationship  Among  Different  Coordinate  Systems 

2.3  Aircraft  Modeling 

A  model  of  an  aircraft  is  a  mathematical  description  of  the  aircraft  motion.  It  is  given  by  a  set 
of  differential  equations,  called  the  equations  of  motion,  which  present  the  relations  between 
the  control  inputs  and  the  state  variables.  From  the  equations  of  motion,  we  hope  to  gain 
some  physical  intuitions  that  will  help  to  structure  the  model  in  a  real-time  simulator. 

The  equations  of  motion  are  derived  based  on  first  principles.  Detailed  derivations  and  vari¬ 
able  definitions  are  given  in  Appendix  B.  Since  landing  is  the  main  concern  in  this  report,  we 
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will  consider  the  flat-earth  model  derived  in  Appendix  B  and  rewrite  the  equations  of  motion 


as  follows: 

Pned  = 

R>B 

F 

-£lB\B  +  Ran&o 

m 

—  J  1  £2  B  Jco  B  +  J  1 

<I> 

where  g0  is  the  gravity  expressed  in  the  NED  frame,  Ffi  andTfi  are  the  resultant  aerody¬ 
namic  and  thrust  forces  and  the  moments  they  generate  expressed  in  the  ABC  frame,  and 

Pned  :  the  position  of  aircraft  eg  expressed  in  the  NED  frame, 

\  b  :  relative  velocity  of  aircraft  eg  w.r.t  air  mass  expressed  in  ABC  frame, 

co  B  :  absolute  angular  velocity  of  ABC  frame  expressed  in  ABC  frame, 

O :  Euler  angle  between  ABC  and  NED  frames. 


For  the  sake  of  linearization,  we  modified  the  above  equations  by  expressing  the  relative  ve¬ 
locity  of  aircraft  eg  and  the  absolute  angular  velocity  of  the  ABC  frame  in  the  ABC- Wind 
frame.  Since  vK,  =  RWBvB,<aw  =  RWBe>B  andFM.  =  RWBFB ,  we  have  the  revised  equations  of 

motion  as 


Navigation  Equations : 
Force  Equations : 
Kinematics  Equations : 
Moment  Equations : 


Pned  ~ranrwbvw 

v„.  -^u  vu  +  rwbran£o  +K  tm 

O  =  ¥(*)/&  <oH. 

= -sc®  »■  -  (*>„■ x  j;X 


(i) 


where 


0 

~P 

-dr  cos  p 

0 

—  0t)Z 
w 

K  ' 

ft 

0 

dr  sin  P 

K 

0 

-K 

dr  cos  P 

~a  sin  (5 

0 

-< 

cox 

w 

0 

Noticing  that 

\w  =  [vr  ,0,0]T  with  vT  the  true  air  speed 
and 

vM,+Q^vlv=[vr,  fivT,  avj  cos  p]T  > 
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we  can  write  the  equations  of  motion  with  the  state  variables  as  follows: 


Pe >  P  >  R  »  2> 

by  letting  pNED  =[Pa,,  p£,/i]r  ,  vlt,  =  [v7.,0,0]r ,  <J>  =  [0,0,p]r ,  and  to =  [P, 2, i?]T  .  Fur¬ 
thermore,  we  specify  the  force  and  torque  in  terms  of  aerodynamic  and  thrust  forces  and  the 
moments  they  generate  as 


'-  D 

'Ft' 

T 

Y 

+  rwb 

0 

1 

ft* 

II 

M 

-L 

0 

N 

L  =  qSbCf  -  Rolling  moment 
M  =  qScCM  -  Pitching  moment 
N  =  qSbCN  -  Yawing  moment 


where  Fr  is  the  thrust, 

D  =  qSCD  -  Drag 
L  =  qSCL  —  Lift 
Y  =  qSCy  -  Side  force 

and 

CD  -  Drag  coefficient 
C L  -  Lift  coefficient 
CY  -Side  -  force  coefficient 
q  -  free  stream  dynamic  presure 
c  -  wing  mean  geometric  chord 


C,  -  Rolling  moment  coefficient 
CM  —  Pitching  moment  coefficient 
CN  -  Yawing  moment  coefficient 
S  -  wing  reference  area 
b  -  wing  span 


The  forces  and  moments  in  Equation  2  are  the  variables  that  affect  the  dynamics  of  the  sys¬ 
tem,  but  they  are  not  the  control  variables.  In  fact,  the  force  and  moment  coefficients  C.,  with 
*  being  the  subscripts  defined  in  Equation  2,  are  functions  of  the  control  variables: 

u  =  [thl,  el,  ail,  rdr]T  (3) 


where 

fW-Throttle  ai/-Aileron  deflection 

^/-Elevator  deflection  rdr-Rudder  deflection 
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Equations  1-3  compose  the  nonlinear  mode  of  the  aircraft  considered  in  this  report.  The  air¬ 
craft  dynamics  are  often  defined  by  the  last  three  sets  of  equations  in  Equation  1  (i.e.,  the 
force  equations,  the  kinematic  equations,  and  the  moment  equations). 

The  aircraft  motion  can  be  further  characterized  as  a  combination  of  longitudinal  motion  and 
lateral-directional  motion.  While  the  former  describes  the  pitching  and  translational  part  of 
the  motion  with  variables  a,0,vT,Q,  the  latter  is  about  rolling,  sideslipping,  and  yawing 
with  variables  /?,  (, p ,  q>,  P,  R .  When  an  aircraft  flies  level,  <p  =  P  =  0 ,  and  non-sideslipping 
(coordinated  flight),  /5  =  R  -  0 ,  the  longitudinal  motion  can  be  described  by  a  set  of  decou¬ 
pled  equations 


mvT  =  FT  cos  or  -  D  -  mg  0  sin  {9 -a) 

mavj  =  -Ft  sin  a  -  L  +  mvTcol  +mg0  cos (6  -  or) 

e=col 

Jy6>i=M 


Linearization  of  the  nonlinear  model  of  an  aircraft  can  be  derived  with  respect  to  a  flight 
condition  (i.e.,  a  description  of  a  particular  flight  situation).  A  commonly  considered  flight 
condition  is  level,  coordinated  flight,  which  is  specified  by 

P,  0,  P,Q,R  =  0,  x  =  0,  u  =  0 


This  flight  condition  implies  the  following  steady  state: 

xs  =[vr,  or,  6,  <f>,0, y/,P , <2i^]l/j,^,p,g,R,x=o=[v7'J >  0,  6S,  y/s,  0,  0,  0] 

(4) 

us  =[thl,el,ail,rdr]\/}4PQRis0=[thls,  els,  ails ,  rdrs] 

Define  Sx  =  x-xs,  Su  =  u-  us ,  then  the  linearized  system  at  the  steady  state  in  Equation  4  is 
given  by 

E8x-  ASx  +  BSu 


where  E,  A,  and  B  are  constant  matrices  which  can  be  found  in  Section  2.5  of  the  book  by 
Stevens  [Stevens  92], 

The  objective  of  analytically  deriving  the  linearized  model  is  to  understand  the  structure  of 
the  aircraft  dynamics.  In  particular,  we  would  like  to  see  what  elements  in  matrices  E,  A,  and 
B  are  intrinsically  zero  or  invariant  as  the  type  of  aircraft  or  the  flight  condition  change. 
Those  elements  will  be  treated  as  the  fixed  parameters  in  system  identification.  As  can  be 
seen  from  above,  analytically  derived  linearization  could  be  tedious,  as  the  process  has  to  be 
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repeated  when  the  flight  conditions  are  changed.  This  motivates  the  development  of  numeri¬ 
cal  linearization,  which  is  discussed  in  Section  2.4. 


2.4  Modeling  Tools 

Stevens  describes  how  insight  into  the  dynamics  of  an  F-16  fighter  was  gained  by  experi¬ 
menting  with  a  standard  model  described,  and  the  simulation  experiment  runs  with  the  main 
Fortran  code  provided  by  Stevens  [Stevens  92],  Two  different  sets  of  programs  were  used, 
one  to  find  a  specific  flying  condition  (and  then  the  steady  state)  and  the  other  to  do  a  nonlin¬ 
ear  simulation  of  landing.  Appendix  C  contains  a  list  of  the  files  and  subroutines  for  the 
original  and  modified  source  codes. 

The  core  of  both  programs  is  a  subroutine  that  contains  a  nonlinear  model  of  the  plane.  To 
find  the  steady  state,  an  optimization  routine  is  needed.  The  simulation  program  uses  an  ordi¬ 
nary  differential  equation  solver  (Runge-Kutta  method). 

The  aircraft  system  is  described  by 
x=/(x,u) 

x  =  [  Vt,  a,  p,  <(>,  e,  v,  p,  q,  r,  x,  y,  z,  pow  ] 
u  =  [  thl,  el,  ail,  rdr] 

Table  1  displays  all  the  state  variables  and  inputs  with  the  corresponding  units  as  they  are 
used  in  the  source  code.  This  is  also  the  order  in  which  they  appeared  in  the  subroutine. 


Vt  (ft/s) 

True  air  speed 

a  (rad) 

Angle  of  attack 

P  (rad) 

Side-slip  angle 

<)>,  0,  vj /  (rad) 

Roll,  pitch,  and  yaw  angles 

p,q,r  (rad/s) 

Roll,  pitch,  and  yaw  angle  rates 

x,y,z  (ft) 

XYZ  coordinates  in  flat-earth  system  frame 

el  (deg) 

Elevator  position  (actuator) 

ail  (rad) 

Aileron  position  (actuator) 

rdr  (rad) 

Rudder  position  (actuator) 

thl  (unit) 

Throttle  position  (actuator) 

pow  (unit) 

Power/thrust 

Table  1:  Aircraft  Nomenclature 
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The  physical  model  must  be  enhanced  with  actuator  dynamics  and  constraints  and  the  use  of 
the  typical  inputs  that  the  pilot  has  available  (i.e.,  thl,  the  roll  command  input  Pcom,  and  the 
pitch  command  Qcom).  This  input  transformation  translates  in  an  augmented  system 

x  =  /(x,u) 

x  =  [  Vt,  a,  (3,  <j),  0,  \\f,  p,  q,  r,  x,  y,  z,  pow,  this,  els,  ails,  rdrs  ] 
u  =  [  thl,  Pcom,  Qcom], 

with  constraints  on  the  actuator  states  this,  els,  ails,  and  rdrs  and  their  derivatives.  It  should 
be  noticed  that  pitch  and  roll  command  inputs  are  actually  indirect  input  commands.  The  air¬ 
craft  is  directly  controlled  by  the  position  of  the  elevator,  ailerons,  and  rudder  flaps.  There¬ 
fore,  there  is  an  implicit  internal  conversion  system,  perhaps  a  feedback  controller,  that 
translates  pitch  and  roll  commands  into  actuator  commands.  In  Section  3,  we  show  an  exam¬ 
ple  of  this  program. 
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3  Automatic  Landing  of  an  F-1 6  Aircraft 


In  this  section,  we  are  concerned  with  the  last  stages  of  the  landing  procedure  (i.e.,  when  the 
aircraft  is  in  relatively  good  alignment  with  the  runway).  The  automatic  landing  system 
should  provide  feedback  control  to  allow  the  aircraft  to  follow  a  trajectory  in  the  XYZ  space. 

The  trajectory,  called  glide-path ,  is  determined  by  the  position  and  orientation  of  the  runway 
and  the  flight-path  angle  yR.  A  typical  automatic  landing  system  uses  a  radio  beam  directed 
upward  from  the  ground  that  can  be  detected  by  aircraft  sensors  to  control  the  landing.  Figure 
4  shows  a  landing  approach  and  some  important  characteristics. 

The  control  commands  to  be  used  are  the  usual  three:  throttle,  pitch,  and  roll  command  in¬ 
puts.  The  throttle  command  is  generally  left  open-loop  and  it  is  directly  controlled  by  the  pi¬ 
lot  or  by  a  separate  subsystem.  Pitch  and  roll  commands  usually  have  some  level  of  feedback 
control  to  help  the  pilot.  The  pilot  can  modify  the  command  actions  using  the  stick. 


Aircraft  trajectory 
\ 


\ 

N 

\ 

Glide-path  ***\ 


s. 


Abort-landing 

^ - 


Decision-height 


N. 


Runaway 


~  10  secs 


Figure  4:  Landing  Trajectoy 

Landing  specifications  include  constraints  on  the  way  the  approach  should  be  taken  (i.e.,  a 
desired  range  of  values  for  the  aircraft  state  variables  and  flight-path).  Table  2  shows  typical 
values  for  a  landing  approach. 
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Flight-path  angle 

constant  between  2.5  to  3.5  degrees 

Airspeed 

240-280  ft/s 

Angle  of  attack 

between  10-15  degrees 

Table  2:  General  Landing  Specifications 

One  of  the  most  important  features  of  an  automatic  landing  system  involves  the  decision  of 
whether  to  land  the  aircraft  according  to  how  the  approach  has  been  going  up  to  that  time.  A 
way  to  address  this  issue  is  by  defining  a  position,  taking  into  account  the  aircraft  dynamics, 
when  the  decision  to  land  or  abort  has  to  be  made.  This  point  is  called  decision-height. 
Tighter  constraints  on  the  state  of  the  aircraft  need  to  be  met  at  decision-height  to  proceed 
with  the  landing.  Table  3  specifies  a  typical  set  of  values  that  are  required  at  the  decision- 
height  point  to  continue  and  land  for  the  state  of  an  F-16  aircraft.  Actually,  the  specification 
requires  that  these  constraints  be  satisfied  from  the  decision-height  to  the  touchdown  point. 


Decision-height  point 

roughly  10  seconds  before  landing 

Angle  of  attack  (a) 

between  10-15  degrees 

Vertical  deviation  (dv) 

maximum  5  feet 

Horizontal  deviation  (dh) 

maximum  15  feet 

Roll  error  angle  (<()) 

Maximum  5  degrees 

Pitch  error  angle  (0) 

Maximum  5  degrees 

Yaw  error  angle  (\| /) 

Maximum  5  degrees 

Sink  rate 

Between  250  and  1000  ft/min 

Table  3:  Typical  Range  of  Values  for  Aircraft  Variables  at  Decision-Height  Point 

Rather  than  using  the  absolute  position  of  the  aircraft  in  the  XYZ  coordinate  system,  the  error 
distance  to  the  glide  path  is  used  to  evaluate  the  approach  and  design  the  control  system.  Two 
scalar  values  are  considered:  the  vertical  and  horizontal  deviations  dv  and  dh,  and  the  vertical 
and  horizontal  projected  distance  of  the  aircraft  with  respect  to  the  glide  path.  Angular  mag¬ 
nitudes  can  be  considered  instead,  but  dv  and  dh  give  easier  equations  with  which  to  work. 
Having  set  up  this  goal  with  respect  to  the  decision-height  point,  the  next  step  is  to  design  a 
controller  that  brings  the  aircraft  into  that  region.  A  similar  goal  is  to  assess  the  performance 
of  a  previous  controller  and  predict  when  or  where  that  controller  would  take  us  to  a  good 
landing.  These  problems  are  addressed  in  the  next  two  sections. 
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4  Model  Identification 


The  objective  of  this  section  is  to  identify  a  model  for  the  F-16  simulator  using  the  data  col¬ 
lected  while  it  is  simulating  landing.  Rather  than  a  black-box  identification,  it  proves  very 
useful  to  use  the  information  from  a  first-principles  model  and  standard  control  techniques  to 
establish  the  best  structure  for  the  model  to  be  identified.  The  overall  aircraft  control  system 
can  be  illustrated  as  in  Figure  5,  where  thl,  Pcom,  and  Qcom  are  the  control  variables  for 
which  the  control  algorithms  will  be  designed.  Together  with  the  inner  loop  control,  these 
control  variables  generate  control  commands  for  the  physical  control  surface  thl,  el,  ail,  and 
rdr,  which  in  turn  control  the  aircraft.  The  model  that  we  need  to  identify  consists  of  the  air¬ 
craft  dynamics,  actuator  dynamics  and  inner  loop  control.  We  will  proceed  in  two  steps.  As 
studied  in  Stevens,  the  first  step  is  to  get  the  F-16  model  by  running  the  software  provided  by 
Stevens,  and  then  the  second  step  is  to  use  the  structure  of  the  model  obtained  for  the  simu¬ 
lator  and  identify  the  unknown  parameters  [Stevens  92]. 


thl 

Internal 

— I — ^ 

W 

Pcom 

Controller 

- ► 

Actuator 

w 

Qcom 
- w. 

(inner  loop 
control) 

- ► 

Dynamics 

W 

- ^ 

I 


thl 

~ef 

ail 

rdr 


> 

>  Aircraft 
^  Dynamics 

>  _ 


x 


Figure  5:  F-16  Control  System 

4.1  A  Book  Model  for  F-16  Aircraft 

We  used  the  nonlinear  model  studied  in  the  book  by  Stevens  to  analyze  the  characteristics  of 
the  aircraft  system  and  propose  a  sound  structure  for  the  linear  model  [Stevens  92].  Since  we 
propose  a  control  design  based  on  a  linear  model,  we  need  to  find  a  linearized  model  for  the 
aircraft  in  landing  conditions.  This  was  done  by  using  the  Fortran  programs  provided  by 
Stevens,  specifically  the  trim  program.  After  that,  we  analyzed  possible  structures  for  the  in¬ 
ner  loop  controller  that  meets  the  demands  of  the  landing  specifications.  Simulations  were 
carried  out  using  a  modified  simu  program,  the  book  simulator. 

4.1.1  Trim  Data 

Based  on  the  landing  specifications  discussed  in  Section  2,  the  If  16  trimming  program  was 
used  to  generate  a  steady  state  and  linearize  the  six-degrees  of  freedom  nonlinear  model  em- 
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bedded  in  the  program.  Figure  6  and  Figure  7  show  a  listing  of  the  input  and  output  of  the 
trimming  program,  respectively. 


DRIVER  PROGRAM  FOR  TRIMMING  &  LINEARIZING  Fl6 
eg  position  ?  (def.=  0.B5)  :  0.3 

Trim  (T)  or  Linearize  (L)  :  t 

****  SUBROUTINE  TO  FIND  STEADY-STATE  FLIGHT  CO 
Copyright  Jan.  1992,  B.  L.  Stevens 


Enter  Aircraft  States  : 

Position  coords  (ft):  North,  East,  &  Alt.:  0,0,0 

Climb  angle  &  Compass  Heading  (deg):  -2.5,0 

True  air  speed  (ft/s):  260 

Roll,  Pull-up,  &  Turn  rates  (degs/s) :  0,0,0 

Reqd.  #  of  trim  iterations  (def.  =  1000):  10000 


Figure  6:  Input  Listing  for  the  Trimming  Program 


****  SUBROUTINE  TO  FIND  STEADY-STATE  FLIGHT  CONDITIONS 
Copyright  Jan.  1992,  B.  L.  Stevens 


Enter  Aircraft  States: 

Position  Coords  (ft):  North,  East,  &  Alt.:  0,0,0 

Climb  angle  &  Compass  Heading  (deg):  -2.5,0 

True  air  speed  (ft/s):  260 

Roll,  Pull-up,  &  Turn  rates  (degs/s):  0,0,0 

Reqd.  #  of  trim  iterations  (def.  =  1000):  10000 

Throttle  Elevator,  Ailerons,  Rudder 
1.01E-01  -4.03E+00  -9.56E-07  1.81E-06 


Angle  of  attack  1.21E+01  Sideslip  angle  8.32E-07 

Pitch  angle  9.64E+00  Bank  angle  0.00E+00 

Normal  acceleration  9.86E-01  Lateral  acceln  -1.47E-08 
Dynamic  press.  8.03E+01  Mach  2.33E-01 


initial  cost  function  1.92E+02  Final  cost  function  1.35E-1 


Figure  7:  Output  of  the  Trimming  Program  to  Screen 

Saving  these  conditions  in  a  file,  the  same  program  generates  a  linearization  of  the  nonlinear 
model  around  these  conditions.  The  linearized  model  is  given  in  the  following  form: 


x  =  Ax  +  Bu 
y  =  Cx  +  Du 


The  matrices  A,  B,  C,  and  D  generated  here  allow  us  to  evaluate  which  interactions  are  im¬ 
portant  for  establishing  the  structure  of  the  model  to  match  with  the  model  in  the  real-time 
simulator.  Figures  8  and  9  show  an  example  of  the  output  of  the  trim  and  linearization  results 
with  some  modifications  to  make  them  easier  to  understand.  It  is  important  to  note  at  this 
point  the  strong  structure  of  the  matrices  A  and  B  that  come  out  of  the  linearization  program. 
They  have  many  zeros  establishing  the  decoupling  between  lateral  and  longitudinal  motions 
for  this  particular  case.  Different  conditions  were  used  for  the  center  of  gravity,  climb  angle, 
and  air  velocity  with  similar  results. 
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Figure  8:  Output  of  Trimming  Program  to  a  File 
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Figure  9:  Example  Output  from  Linearization  Program  Showing  Matrices  A  and  B 

From  this  analysis  we  concluded  that  there  was  no  reason  to  try  to  estimate  all  of  the  ele¬ 
ments  in  the  matrices  that  define  the  linear  approximation  model.  Moreover,  we  determined 
which  coefficients  dominated  the  dynamic  behavior  and  needed  to  be  estimated  using  the 
experimental  data. 


4.1.2  Inner  Loop  Control  Structure 

In  order  to  progress  in  the  study  and  get  a  better  understanding  of  aircraft  dynamics,  an  inner 
loop  controller  as  shown  in  Figure  5  had  to  be  designed  for  the  book  simulator.  Classical  and 
modem  control  techniques  could  be  used  for  such  a  task  [Stevens  92].  However,  we  preferred 
to  use  modem  control  techniques  that  allowed  us  to  introduce  the  constraints  and  a  priori 
limits  of  the  linear  model  into  the  picture. 

For  the  landing  conditions,  it  was  fairly  clear  from  the  trim  example  that  lateral  and  longitu¬ 
dinal  dynamics  could  be  separated  in  design  of  the  internal  controller  and  leave  the  possible 
coupling  terms  for  the  higher  levels  of  control.  Also,  the  inner  loop  controller  does  not  have 
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to  deal  with  any  specifics  of  the  landing  route,  so  the  deviations  from  the  glide  path  were  not 
considered  for  feedback.  However,  it  was  necessary  to  introduce  feedback  from  the  angle 
positions  to  control  the  aircraft  with  a  reasonable  behavior. 

A  different  problem  was  posed  by  the  throttle  command.  In  the  Fortran  simulator,  we  had 
complete  control  of  the  aircraft  and  could  design  a  throttle  control  in  the  same  way  that  we 
designed  the  directional  feedback  control.  However,  in  the  real-time  simulator  and  because  of 
flight  specifications,  we  did  not  have  control  over  the  throttle,  which  introduced  another  dis¬ 
turbance  and  internal  dynamics  into  the  control  design.  Therefore,  to  make  things  similar  in 
both  environments,  we  designed  a  decoupled  controller  for  the  throttle  and  worked  with  that 
controller  throughout  the  rest  of  the  study. 

Thus,  the  structure  of  the  inner  loop  controller  chosen  is  as  follows: 

the = thsl  +kth  (Vt-Vts, ) 

elc=elsl  +kel{e-6st)+ke2  (q-qsl  )+ke3  Qcom 

alc=alsl  +ka]{</>-</>s,)+ka2{y/-yfsl)+kai(p-pst  )+ka4(r-rsl  )  +  ka5  Pcom 
rdc=rdsl+kr]{<t>-<psl)+kr2{y/-i//sl)+kri(p-ps,)  +  kr4{r-rst)  +  kr5  Pcom 

with  the,  elc,  ale,  and  rdc  being  the  respective  throttle,  elevator,  aileron,  and  rudder  control 
commands  that  feed  the  actuator  dynamics.  This  controller  structure  was  verified  during  the 
minimization  procedure  involving  LMIs  and  is  explained  in  the  next  section. 


4.1.3  Outer  Loop  Controller 

To  be  able  to  meet  the  specifications  of  the  automatic  landing  procedure,  an  outer  loop  con¬ 
troller  was  designed.  It  introduced  feedback  from  the  glide  path  deviations  to  make  the  neces¬ 
sary  corrections.  Again,  the  minimization  procedure  semidefinite  programming  (SDP) 
showed,  as  expected,  that  the  vertical  deviation  error  is  important  only  for  the  elevator  con¬ 
trol  command  and  that  the  horizontal  deviation  is  significant  only  for  the  aileron  and  rudder 
commands.  The  outer  loop  controller  structure  is  of  the  following  form: 

Qcom  =  kql{0-6st)+kq2  (q-qsl  )  +  kq3  dv 

Pcom  =  kp]{<p-<t>sl)+kp2(\ff-iys,)+kp3  (p- pst  )+kp4  (r-rst  )  +  kp5  dh 

4.2  A  Linear  Model  of  the  Real-time  F-16  Simulator 

In  this  subsection,  we  apply  the  results  obtained  in  the  last  subsection  to  identify  a  linear 
model  for  the  real-time  F-16  simulator  in  landing.  In  particular,  we  use  the  obtained  aircraft 
structure,  the  actuator  characteristics,  and  the  inner  loop  control  to  build  an  augmented  sys¬ 
tem  for  the  real-time  simulator  with  unknown  parameters.  Then  identification  of  the  system 
becomes  an  issue  of  parameter  identification. 
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4.2.1  Model  Structure 

According  to  the  results  from  Section  4.1.3,  we  construct  the  model  as 

x  =  Ax  +  Bu,  y  =  Cx 

with  x  =  [8</>,  89,  S<p,  P,  Q,  R,  8el ,  8ail ,  Srdr  ] , 
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where  S  *  or  8,  are  the  differences  between  the  actual  variables  and  their  steady  state  value, 
and  parameters  a  -  and  by  in  A  and  B  matrices  are  to  be  identified. 

The  physical  interpretation  of  this  structure  is  clear.  In  fact,  partitioning  the  state  x  as 
x  =  [xx,x2]  with  =[<50,  89,  8(f>,  P,  Q,  R]  and  x2  =  [8eI ,  Sa!l ,  8rdr ] ,  and 


A  = 

"Ai 

A  2 

,  B  = 

% 

o 

t 

_^21 

a22_ 

l*2j 

accordingly,  it  can  be  seen  that  matrices  A} ,  and  An  determine  the  aircraft  dynamics,  A21 
constructs  the  inner  loop  controller,  and  A22  represents  the  actuator  dynamics.  The  block 
diagram  shown  in  Figure  10  illustrates  such  an  interpretation. 
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Figure  10:  A  Block  Diagram  for  the  Structure  of  the  Model 

4.2.2  Parameter  Identification 

To  identify  the  parameters,  we  formulated  a  minimization  problem  and  used  the  Matlab 
function  fmins  to  solve  it.  The  detailed  procedure  is  described  in  the  following  steps. 

Step  1:  Collect  data  from  the  real-time  simulator.  Specifically,  fly  the  aircraft  through  the 
landing  stage  and  record  the  values  of  the  state  variables.  The  state  variables  used  in  pa¬ 
rameter  identification  are 

state  _  data  =  [</>,  9 ,  tp,  P,  Q,  R)  -  steady  _  state  _  state 
and  the  control  variables  are 

control  _  data  =  [Qcom ,  Pcom  ]  -  steady  _  state  _  control 

where  steady _state_state  and  steady _state_control  are  the  mean  values  of  the  involved  vari¬ 
ables.  It  is  important  that  the  data  are  obtained  from  a  successful  landing. 

Step  2:  Solve  a  minimization  problem  using  Matlab  function  fmins.  Tire  function  fmins  is 
used  in  the  form 

par  =fminsC  F_name\  par0 ) 
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which  attempts  to  return  a  vector  par  that  minimizes  objective  function  F(par)  defined  in  the 
file  F_name.m.  par0  is  an  initial  guess  of  the  values  of  par.  To  use  function  fmins  in  parame¬ 
ter  identification,  we  define  the  objective  function  as  the  absolute  error  between  the  data  col¬ 
lected  from  the  simulator  and  the  trajectory  generated  by  the  model  with  trial  parameters. 
Specifically,  for  each  set  of  parameters  (starting  with  an  initial  guess),  we  complete  the  ma¬ 
trices  A  and  B,  simulate  the  trajectory  of  system  x  =  Ax  +  Bu  with  the  initial  condition 

x(0)  =  [0(ro ),  0(to ),  (p{t0  ),P(t0  ),Q(t0  )JR(t0 ),  0, 0, 0]  -  [steady  _  state  _  state,  0, 0, 0] , 

and  control  input  control_data,  and  then  compute  the  error 

error  =*x(l)-</>\  +  \x(2)-0\  +  \x(3)-<p\  +  \x(4)-P\  +  \x(.5)-Q\  +  \x(6)-R\. 

After  the  execution  of  fmins  terminates,  vector  par  contains  the  parameters  to  be  identified. 

Step  3:  Verify  the  obtained  model  by  comparing  the  model  trajectories  with  the  real  data  ob¬ 
tained  from  another  run  of  the  landing  simulation.  If  the  result  is  satisfactory,  model  identifi¬ 
cation  is  complete.  Otherwise,  try  another  set  of  initial  parameters  and  repeat  Step  2. 

We  now  apply  the  described  procedure  to  parameter  identification.  As  mentioned  earlier,  the 
model  that  we  will  identify  is  with  respect  to  a  particular  flight  condition  (i.e.,  landing  in  this 
case).  Hence,  it  is  important  to  make  sure  that  the  real  data  we  use  are  indeed  collected  when 
the  aircraft  is  landing.  Figure  1 1  shows  the  altitude  of  the  aircraft  in  the  group  of  data  gath¬ 
ered. 


Figure  11:  Altitude  of  the  Aircraft  in  Landing 

As  can  be  seen  from  Figure  11,  the  aircraft  was  actually  landing  after  120  seconds.  Therefore, 
the  portion  of  data  to  be  used  for  identification  should  be  between  time  120  and  190  seconds. 
With  this  selected  set  of  data,  the  parameters  are  identified  by  working  through  Step  2  in  the 
procedure.  As  a  result,  we  obtain  the  model  given  as  follows: 
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Figure  12  shows  the  trajectory  of  the  identified  system  plotted  against  the  real  data  used  for 
identification  in  the  time  interval  [145,  190]  seconds.  It  can  be  seen  that  the  simulated  data 
match  the  real  data  reasonably  well.  Figures  13  and  14  illustrate  the  model  verification  result, 
where  Figure  13  depicts  the  result  for  the  landing  period  and  Figure  14  displays  the  compari¬ 
son  for  the  whole  range  of  data.  Again,  the  match  between  the  simulated  data  and  the  real 
data  is  good  except  for  some  disturbance  in  the  real  data.  In  particular,  the  result  in  the  long 
run  demonstrates  that  the  model  identified  is  reasonably  robust  (i.e.,  it  represents  the  real 
system  fairly  well  even  with  variations  of  flight  condition). 


During  landing,  the  aircraft  should  be  controlled  so  it  flies  along  the  glide  path.  Hence  the 
control  problem  becomes  stabilization  of  the  aircraft  on  the  glide  path.  In  the  set  of  measured 
data,  the  deviations  from  the  glide  is  described  by  two  variables,  gv  and  gh,  the  angles  be¬ 
tween  the  straight  line  from  aircraft  eg  to  the  touch-down  point  and  the  glide  path  in  the  ver¬ 
tical  plane  and  horizontal  plane,  respectively.  Figure  15  shows  these  two  angles. 


Roll  angle  PHI  Pitch  angle  THETA  Yaw  angle  YAW 


Roll  rate  P 


Pitch  rate  O 


Yaw  rate  R 


Figure  12:  Trajectory  of  Identified  Model  (Solid  Line)  and  Real  Data  Used  for 
Identification  (Dotted  Line) 
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Figure  13:  Results  of  Model  Verification  in  the  Landing  Period  (Solid  Line:Trajectory 
of  the  Identified  Model,  Dotted  Line:  Real  Data) 
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Figure  14:  Results  of  Model  Verification  in  a  Long  Run  (Solid  Line:  Trajectory  of 
Identified  Model,  Dotted  Line:  Real  Data) 
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Figure  15:  Angles  Determining  Glide  Path  Derivation 

When  a  control  algorithm  is  designed  to  stabilize  the  aircraft  along  the  glide  path,  it  is  prefer¬ 
able  to  regulate  the  aircraft  deviation  from  the  glide  path,  namely  the  distance  from  the  air¬ 
craft  eg  to  the  glide  path.  Such  a  distance  can  be  resolved  into  two  components,  the  horizon¬ 
tal  deviation  dh  in  the  horizontal  plane  parallel  to  the  Earth’s  surface,  and  the  vertical 
deviation  dv  in  the  vertical  plane  perpendicular  to  the  Earth’s  surface  and  containing  the 
glide  path.  Then  referring  to  Figure  16,  we  obtain  the  vertical  deviation  to  the  glide  path  as 
follows: 


dv  =VT  sin(/-  yR), 


d=- 


sm  gv 


sin  i-yR+ g  v) 


and  the  horizontal  derivation  to  be 
dh  =VTsin(<p-<pR),  dh  =<Jpl+pl  sing*, 
where  VT  is  the  steady-state  value  of  the  air  speed  vT  . 


(5) 


(6) 
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Figure  16:  Derivations  to  the  Glide  Path 

To  incorporate  the  derivations  to  the  glide  path  in  the  model,  we  linearize  the  differential 
equations  of  in  Equations  5  and  6  as 

dv=VTSe,  dh=VTS<p. 

Here  we  have  assumed  that  the  air  speed  vT  and  the  angle  of  attack  a  are  kept  constant  by 
the  throttle  control,  and  the  steady-state  value  vT  is  VT  =  260  ft/sec .  Then  the  complete 
model  to  be  used  in  control  design  involve  11  states,  which  are  as  follows: 

x  =  [S</>,  56,  5<p,  P,  Q,  R,  dv,dh,  5el,  Sail ,  5rdr ) 
and  matrices  A  and  B  are  given  by  the  following: 
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This  will  be  the  model  considered  in  the  next  section  when  the  control  algorithm  is  designed. 
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5  Stabilization  Control  and  Stability 
Analysis 


Based  on  the  linearized  model  developed,  control  algorithms  were  designed  to  stabilize  the 
aircraft  along  the  glide  path.  Since  the  designed  control  laws  were  implemented  as  the  base¬ 
line  control,  we  also  derived  the  stability  region  of  the  closed-loop  system  as  the  criterion  for 
switching  to  the  baseline  controller  when  the  upgraded  controller  causes  the  aircraft  to  be¬ 
have  abnormally.  In  this  section,  we  review  the  LMI-based  control  design  and  stability  deri¬ 
vation  discussed  in  Seto’s  technical  report,  then  apply  the  approach  to  the  stabilization  prob¬ 
lem  for  F-16  landing  [Seto  99].  Finally,  we  present  some  experiments  with  the  book 
simulator,  where  the  proposed  design  approach  is  used  to  design  various  control  algorithms 
for  the  book  model  of  the  F-16  in  [Stevens  92]. 

5.1  LMI-Based  Control  Design  Approaches 

In  this  subsection,  LMI  techniques  are  used  to  design  the  feedback  controller  with  state  and 
input  constraints  [Boyd  94].  Since  the  performance  of  the  controlled  system  is  determined  by 
the  poles  of  the  closed-loop  system,  constraints  on  the  pole  location  in  the  complex  plane 
may  be  specified  to  enforce  certain  performance.  Such  constraints  can  also  be  expressed  in 
terms  of  LMI  [Chilali  96].  Details  are  provided  below.  Consider  a  linear  system  in  state  space 
form 

x  =  Ax  +  Bu,  xeRn,ueRm 

subject  to  the  state,  control,  and  rate  constraints  as 

djX<\,  a^R",  i  =  l,...,q,  bjU  <  1,  bj€Rm,  j  =  l,...,r,  and  ckx<l,k  =  1, 

Our  main  goal  was  to  design  a  state  feedback  controller  u  =  Kx  that  meets  the  constraint 
specifications  and  find  out  its  maximal  ellipsoidal  stability  region 

S  =  {x:  xTPx<  1} 

with  P  a  symmetric  positive  definite  matrix,  denoted  by  P  >  0. 
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Stability,  constraints,  and  pole  placements  can  all  be  expressed  in  terms  of  LMI  in  this  case. 
Stability  reduces  the  existence  of  an  nx  n  symmetric  matrix  Q  and  an  m  x  n  matrix  Y  that 
satisfy  the  following  LMI  constraints: 

0>  0 

(8) 

AQ  +  QAT  +BY  +  YtBt  <0 

Then  the  control  gain  is  given  by  K  =  YQ~' ,  and  the  ellipsoid  is  determined  by  P  =  Q~' .  The 
constraints  on  the  states,  at  x  <  1 ,  are  translated  into  the  following  LMIs: 

diQaJ  <1,  i  =  \,...,q. 

The  constraints  on  the  controls,  bjU  <  1 ,  can  be  expressed  in  LMI  as 


1  bjY 

YTb ]  Q 


>0,  j  =  l,...,r , 


and  the  constraints  on  the  rate  of  change  of  the  states,  ckx<  1 ,  can  be  written  in  LMI  as 


1  ckAQ  +  ckBY 


(ckAQ  +  ckBY) 


>0,  *  = 


The  location  of  the  poles  of  the  closed-loop  system  can  be  restricted  to  a  prescribed  region  in 
the  complex  plane,  such  as  semi-planes,  disks,  or  cones,  to  meet  the  performance  require¬ 
ment.  For  example,  to  reach  a  decay  rate  greater  than  a  in  the  closed-loop  response,  the  poles 
have  to  be  in  the  complex  semi-plane  R e(pole)  <  -a .  This  is  easily  realized  by  modifying 

Equation  8  as  follows: 

<2>0 

AQ  +  Q A'  +BY  +  Y’  B'  +2aQ<0 


For  the  poles  to  be  inside  a  disk  of  center  -q  and  radius  r,  the  following  LMI  has  to  be  added 
to  Equation  8: 


-rQ  qQ+AQ 
qQ+QA'  -rQ 


When  the  desired  region  is  a  cone  with  its  center  in  the  origin  of  the  complex  plane  and 
symmetric  to  the  real  axis  with  angle  0,  the  following  LMI  must  be  specified  in  addition  to 
Equation  8: 
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"sin 6{AQ+QAt)  cos 0(AQ-QAt)  <q 

_cos6(QAt  -AQ  )  sm0(AQ  +  QAT  )J 

These  constraints  offer  good  options  to  set  up  the  pole  locations  for  the  closed-loop  system 
because  it  allows  us  to  set  up  a  decay  rate  and  to  limit  overshoot  and  settling  times.  Figure  17 
displays  a  possible  combination  of  the  LMI  constraints  described  above. 


Figure  1 7:  Pole  Placement  Constraint  Options 

To  maximize  the  ellipsoid,  we  have  to  maximize  det(Q )  or  equivalently  minimize  -log 
det( Q).  The  computation  of  the  optimization  problem  is  performed  using  the  SDPSOL  pack¬ 
age  developed  by  Shao-Po  Wu  and  S.  Boyd  [Wu  96].  In  the  subsequent  subsections,  we  for¬ 
mulate  the  stabilization  problem  for  F-16  landing  as  an  optimization  problem  in  terms  of 
LMIs  and  solve  it  with  the  SDPSOL  software  package.  We  first  carry  out  the  procedure  with 
the  identified  model,  and  then  apply  it  to  the  book  model  and  test  the  rest  with  the  book 
simulator.  Pole  location  constraints  are  introduced  to  the  control  design  with  the  book  model. 

5.2  Control  Design  for  the  Identified  F-16  Model 

We  now  apply  the  approach  described  above  to  the  design  of  a  state  feedback  control  such 
that  the  aircraft  will  fly  along  the  glide  path  and  the  stability  region  of  the  closed-loop  system 
is  maximized.  Since  the  aircraft  is  in  the  landing  mode,  there  are  certain  constraints  on  the 
state  of  the  aircraft,  called  flight  envelope  that  must  be  satisfied.  Such  an  envelope  restricts 
the  aircraft’s  dynamics  to  a  greater  extent  when  it  is  approaching  the  touchdown  point.  This 
requires  the  controller  to  be  more  sophisticated  in  term  of  providing  precise  control  along  the 
glide  path  as  the  aircraft  approaches  the  runway.  Therefore,  we  will  design  two  controllers: 
one  will  be  used  before  the  aircraft  reaches  the  decision  height,  and  another  will  be  in  charge 
after  the  aircraft  reaches  the  decision  height.  The  former  will  be  designed  such  that  the  sta¬ 
bility  region  of  the  closed-loop  system  is  larger  than  the  latter  controller’s  stability  region,  but 
the  former  may  take  a  longer  time  than  the  latter  to  drive  the  aircraft  to  converge  to  the  glide 
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path.  These  controllers  are  designed  by  solving  a  LMI  problem  described  earlier  for  each, 
with  different  sets  of  state  constraints.  Consider  a  system  in  the  form 

x  =  Ax  +  Bu 


with  matrices  A  and  B  given  in  Equation  5  and 
x  =  [S<p,  86,  8<p,  P,  Q,  R,  dv,  dh,  8el,  Sai, ,  8rdr],  u  =  [Sqcom,  Spcom] 


subject  to 


lx,  l<20deg 

lx4  l<  40  deg/  sec 

O 

IT) 

V 

r- 

* 

lx9  l<20deg 

lx9 

1  <  60  deg/  sec 

\x2  l<20deg 

lx5  l<  40  deg/  sec 

X 

oo 

A 

O 

1  x10  l<20deg 

,  I<  80  deg/ sec 

|jt3  l<20deg 

lx6 1<  40  deg/  sec 

1  x, ,  l<20deg 

1  xn 

k  120  deg/ sec 

before  the  aircraft  reaches  the  decision  height  (before  DH),  and 

lx,  l<5deg 

1  x4  l<10deg/sec 

1  x7  l<  5  ft 

lx9  l<24deg 

1  Xg 

l<  60  deg/ sec 

1  x2 1  <  2.5  deg 

lx5  l<10deg/sec 

1  xg  1  <  1 5ft 

1  xl0  l<20deg 

1  -^lO 

l<  80  deg/ sec 

lx3  l<5deg 

lx6  l<10deg/sec 

1  x, ,  1  <  29  deg 

1  ^1 

1  <  1 20  deg/ sec 

after  the  aircraft  reaches  the  decision  height  (after  DH).  The  control  objective  is  to  find  the 
control  gain  K  for  the  linear  state  feedback  control  u  =  Kx  such  that  the  closed-loop  system 
x  =  (A  +  BK)x  is  asymptotically  stable  with  the  above  constraints  satisfied  and  its  stability 

region  is  maximized. 


The  solution  to  the  problem  is  obtained  by  solving  the  symmetric  matrix  Q  and  matrix  Y  from 
the  following  LMI  problem: 


minimize 

log  det  Q 

subject  to 

Q>  0; 

AQ  +  QA 7  +BY  +  YtBt  <0; 

ciiQaJ  <1, 

/  =  11 

1 

( ckAQ  +  ckBY)T 


ckAQ  +  ckBY 

Q 


>o, 


it  =  1,2,3 


where  at ,  i  = 1 ,  are  the  rows  of  the  matrix 


diag(l/20,  1/20,  1/20,  1/40,  1/40,  1/40,  1/50,  1/50,  1/24,  1/20,  1/29) 


before  DH  or 
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diag(l/5,  1/2.5,  1/5,  1/10,  1/10,  1/10,  1/5,  1/15,  1/24,  1/20,  1/29)  after  DH, 


and  c,  =[0Ix9, 1,0,0],  c2  =  [0lx9, 0,1,0],  c3  =  [0lx9, 0,0,1] .  Therefore,  the  complete  solution  to 
the  problem  will  be  two  sets  of  control  gains  with  two  corresponding  stability  regions.  Using 
the  SDPSOL  software  package  [Wu  96],  we  get  the  solutions  to  these  two  LMI  problems  as 
follows: 

Before  the  aircraft  reaches  the  decision  height,  the  control  K  should  be  as  follows: 

-0.0000-0.2589  -0.0000  -0.0000  -0.1745  -0.0000  -0.0156  -0.0000  -0.1582  -0.0000  -0.0000 
0.0207-0.0000  -0.4737  0.7265  -0.0000  -0.5042  -0.0000  -0.0217  0.0000  0.3186  0.3200 

and  the  symmetric  matrix  Q  is  as  follows: 


400.0000 

-0.0008 

-17.1723 

-211.8347 

0.0008 

24.4667 

-0.0317 

-119.9860 

-0.0035 

-59.7096 

113.7391 

-0.0008 

111.2293 

0.0002 

-0.0006 

-70.6586 

-0.0002 

-198.2237 

0.0008 

31.5304 

-0.0021 

-0.0012 

-17.1723 

0.0002 

48.7876 

9.7757 

-0.0002 

-21.7717 

-0.0022 

-136.9569 

0.0008 

-8.4224 

-7.4897 

-211.8347 

-0.0006 

9.7757 

436.7589 

-0.0003 

-30.2084 

0.0030 

17.4860 

0.0029 

-20.5461 

-138.9358 

0.0008 

-70.6586 

-0.0002 

-0.0003 

67.8700 

0.0002 

-0.5276 

0.0007 

43.0385 

0.0003 

0.0006 

24.4667 

-0.0002 

-21.7717 

-30.2084 

0.0002 

14.3932 

0.0005 

10.1714 

-0.0005 

-13.1903 

1.5514 

-0.0317  - 

198.2237 

-0.0022 

0.0030 

-0.5276 

0.0005 

2500.0000 

-0.0724 

-84.7781 

0.0187 

-0.0003 

-119.9860 

0.0008 

-136.9569 

17.4860 

0.0007 

10.1714 

-0.0724 

2500.0000 

-0.0014 

105.3920 

14.7288 

-0.0035 

31.5304 

0.0008 

0.0029 

43.0385 

-0.0005 

-84.7781 

-0.0014 

576.0000 

-0.0026 

-0.0014 

-59.7096 

-0.0021 

-8.4224 

-20.5461 

0.0003 

-13.1903 

0.0187 

105.3920 

-0.0026 

400.0000 

135.0018 

113.7391 

-0.0012 

-7.4897 

-138.9358 

0.0006 

1.5514 

-0.0003 

14.7288 

-0.0014 

135.0018 

663.8984 

After  the  aircraft  reaches  the  decision  height,  the  control  gain  K  is  changed  to  the  following: 

0.0000  -1.4162  -0.0001  0.0000  -0.9437  -0.0001  -0.1570  -0.0000  -0.1521  0.0000  -0.0000 

-0.1313  -0.0000  -1.8460  0.5508  0.0000  -2.9336  -0.0000  -0.0921  -0.0000  0.3127  0.2981 

and  the  symmetric  matrix  Q  becomes  the  following: 


25.0000 

0.0001 

-1.3287 

-13.3985 

-0.0001 

1.7656 

0.0009 

-11.7317 

0.0027 

-19.7653 

1.5411 

0.0001 

4.1497 

-0.0000 

-0.0000 

-4.5740 

0.0000 

-3.7298 

0.0001 

7.4630 

-0.0007 

-0.0005 

-1.3287 

-0.0000 

5.0008 

-1.4746 

-0.0000 

-2.1862 

-0.0001 

-14.3404 

0.0002 

-1.4237 

0.1915 

-13.3985 

-0.0000 

-1.4746 

48.7895 

0.0004 

0.0699 

-0.0003 

4.0581 

0.0021 

-25.9298 

-15.1924 
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-0.0001  -4.5740 

-0.0000 

0.0004 

9.1149 

0.0000 

-1.2334 

0.0001 

20.3156 

0.0019 

0.0010 

1.7656  0.0000 

-2.1862 

0.0699 

0.0000 

1.3932 

0.0001 

1.2728 

0.0004 

-4.5244 

-5.9309 

0.0009  -3.7298 

-0.0001 

-0.0003 

-1.2334 

0.0001 

25.0000 

-0.0046 

-13.6189 

-0.0017 

0.0000 

-11.7317  0.0001 

-14.3404 

4.0581 

0.0001 

1.2728 

-0.0046 

225.0000 

-0.0016 

22.5049 

6.1590 

0.0027  7.4630 

0.0002 

0.0021 

20.3156 

0.0004 

-13.6189 

-0.0016 

576.0000 

-0.0045 

-0.0024 

-19.7653  -0.0007 

-1.4237 

-25.9298 

0.0019 

-4.5244 

-0.0017 

22.5049 

-0.0045 

400.0000 

138.6215 

1.5411  -0.0005 

0.1915 

-15.1924 

0.0010 

-5.9309 

0.0000 

6.1590 

-0.0024 

138.6215 

605.2974 

To  illustrate  the  obtained  result,  we  plot  a  series  of  two-dimensional  projections  of  the  stabil¬ 
ity  region  by  zeroing  out  the  other  nine  variables.  Figure  18  shows  the  results  for  the  case 
before  the  aircraft  reaches  the  decision  height  (larger  stability  region),  and  Figure  19  displays 
the  results  for  the  case  after  the  aircraft  reaches  the  decision  height  (smaller  stability  region). 
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Figure  18:  Two-Dimesional  Projections  of  the  Stability  Region  Before  Aircraft 
Reaches  Decision  Height 
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Q  [-10,10)  deg/sec.  P  [-10,10)  dag/sec. 
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Figure  19:  Two-Dimensional  Projections  of  the  Stability  Region  After  Aircraft 
Reaches  Decision  Height 

To  check  if  an  aircraft  state  is  inside  one  of  the  stability  regions  obtained,  we  need  to  evaluate 
the  function  V  =  xT  Px .  If  the  resulting  V<  1,  the  state  is  inside  the  stability  region;  otherwise 
it  is  not.  In  reality,  variables  S<p,  S6,  Sep,  P,  Q  and  R  are  measurable;  dv  and  dh  can  be  cal¬ 
culated  from  Equations  5  and  6;  but  Sel,  Sail  and  8rdr  may  not  be  available  for  computation 
in  each  sampling  period.  If  this  is  the  case,  an  observer  is  needed  to  estimate  these  variables. 
In  our  development,  we  used  the  identify  observer  in  the  following  form: 

x  =  Ax  +  L(y  -  Cx)  +  Bu 

where  x.  is  the  estimated  state  with  the  last  three  elements  being  the  ones  that  we  are  inter¬ 
ested  in,  matrices  A,  B  and  C  are  defined  in  Equation  7;  u  is  the  control  command  used  in  the 
real-time  simulation;  and  the  observer  gain  L  is  given  by 


16.9301 

-0.0460 

0.0844 

-0.0252 

12.5374 

-0.8626 

0.0349 

-0.8516 

13.3853 

-0.9352 

-0.5823 

1.6391 

0.0397 

-0.0765 

-0.0788 

-0.8829 

-0.6246 

1.4325 

-12.3993 

-4.6792 

21.1685 

-97.3066 

-68.1473 

163.0082 

-144.2297 

-105.5620 

245.3321 

0.9212 

-0.0578 

0.0891 

-0.0170 

0.9634 

-0.0481 

0.0793 

0.0139 

1.1010 

20.4046 

-3.2620 

5.6179 

-1.1621 

20.9738 

-0.3720 

1.4397 

-1.7598 

21.2326 

72.1876 

-391.3307 

28.6433 

272.3444 

-248.7508 

713.2703 

220.5524 

-333.9311 

1143.7588 
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Figure  20  shows  the  comparisons  of  the  real  data  8<p,  86,  8(p,  P,  Q,  and  R  and  the  ob¬ 
served  data,  which  demonstrate  that  the  estimates  match  the  real  data  very  well  after  a  short 
period  of  time. 


PHI 


THETA 


PSI 


P  Q  R 


Figure  20:  Observed  State  (Dotted  Lines)  and  Real  Data  (Solid  Lines) 

To  test  the  obtained  controllers  and  the  stability  regions  before  actual  implementation,  we 
checked  if  the  stability  regions  contain  all  the  states  of  the  aircraft  obtained  from  running  the 
real-time  simulator.  Since  the  data  were  collected  in  a  successful  landing,  the  states  during 
the  landing  should  be  contained  in  the  larger  stability  region  first,  and  in  the  smaller  stability 
region  later.  Figure  21  shows  exactly  this  claim;  namely,  after  120  seconds,  the  aircraft 
started  landing  and  the  value  of  V  was  less  than  1.  From  roughly  165  seconds  on,  the  states  of 
the  aircraft  were  inside  the  smaller  stability  region  until  the  aircraft  touched  down.  The  two 
controllers  and  the  corresponding  stability  regions  were  thus  ready  for  implementation. 


(a)  States  inside  the  larger  stability  region.  (b)  States  inside  the  smaller  stability 

region. 


Figure  21:  Verification  of  the  Aircraft  State  Inside  the  Stability  Regions 
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5.3  Experimental  Tests  with  the  Book  Model 

The  control  design  approach  described  in  Section  5.2  was  tested  with  the  book  simulator. 
Specifically,  we  formulated  the  same  optimization  problem  in  terms  of  the  LMI  constraints 
for  the  book  model,  and  solved  for  the  control  gain  and  the  stability  region.  Then  the  results 
were  implemented  in  the  book  simulator,  and  two  types  of  switching  experiments  were  run  to 
verify  the  feasibility  of  the  approach. 

In  one  type  of  experiment,  two  different  outer  loop  controllers  (recall  the  outer  loop  control¬ 
ler  is  an  implementation  of  the  control  laws  for  qcom  and  pcom)  were  designed  to  work  in  se¬ 
quence;  the  first  one  controls  the  aircraft  up  to  the  decision-height  point,  and  the  second  one 
takes  over  the  control  afterward.  The  first  controller  was  derived  with  relaxed  conditions  on 
the  states  and  glide-path  errors  and  a  decay  rate  to  ensure  that  hard  constraints  were  satisfied 
before  the  decision-height  point  was  reached.  The  second  controller,  on  the  other  hand,  was 
designed  with  restricted  state  and  glide-path  constraints.  Figures  22-25  show  a  switching  ex¬ 
periment  of  this  type. 


Theta  (deg) 


Figure  22:  Pitch  Evolution  During 
Switching  Experiment 


GVdev  (It) 


Figure  23:  Vertical  Derivation 

Behavior  in  a  Switching 
Experiment 


CMU/SEI-99-TR-020 


33 


Figure  24:  Stability  Index  of  Controller  Figure  25:  Stability  Index  of  Controller 
for  Decision  Height  Region  for  Larger  Region 


In  the  other  types  of  experiments,  we  showed  the  recovery  from  a  bad  controller.  Specifically, 
three  controllers  were  running,  and  one  of  them,  say,  Controller  2,  contained  a  bug.  Control¬ 
ler  3  was  designed  to  work  in  a  larger  state  space  than  the  decision-height  specifications,  and 
Controller  1  was  fine-tuned  for  the  decision  height.  Both  Controller  1  and  Controller  3  were 
reliable.  A  stability  represented  by  a  Lyapunov  function  was  derived  with  respect  to  Control¬ 
ler  3,  and  Controller  3  guaranteed  the  stability  of  the  system  so  long  as  the  value  of  the  Ly¬ 
apunov  function  evaluated  at  the  current  state  of  the  system  was  less  than  1 .  As  shown  in 
Figures  26-30,  a  bug  was  introduced  in  Controller  2  at  a  time  equal  to  10  seconds,  which 
eventually  made  the  closed-loop  system  go  out  of  the  stability  region.  Controller  3  took  over 
control  when  it  predicted  the  out-of-bound  Lyapunov  function,  and  it  recovered  the  system. 
At  some  later  time,  Controller  1  was  switched  on  and  landed  the  aircraft. 


Figure  26:  Behavior  of  Angle  of  Attack  Figure  27:  Pitch  Angle  Behavior 

in  a  Recovery  Experiment  During  Recovery 

Experiment 
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Figure  30:  Stability  Indices  for  the  Three  Controllers  During  Recovery  Experiment 
and  Switching  Index  Indicating  Which  Controller  Is  Being  Applied 
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6  Summary  and  Conclusion 


This  report  contains  the  analysis  and  design  of  control  algorithms  for  automatic,  safe  landing 
with  constraint  specifications  for  an  F-16  aircraft.  The  control  objective  of  flying  the  aircraft 
along  the  glide  path  is  formulated  as  a  stabilization  problem  with  constraints.  This  stabiliza¬ 
tion  problem  is  then  translated  to  an  optimization  problem  (that  is,  maximizing  the  stability 
region  subject  to  a  set  of  LMI  constraints,  which  characterize  the  stability  criterion,  the  air¬ 
craft  state  and  control  constraints,  and  the  pole  location  specifications).  This  optimization 
problem  is  solved  by  using  the  SDPSOL  software  package  developed  by  Wu  and  Boyd  [Wu 
96].  As  a  result,  a  set  of  control  gains  for  the  linear  state  feedback  control  law  and  a  stability 
region  are  obtained,  and  both  of  them  are  ready  to  be  implemented  in  the  real-time  simulator. 

While  the  proposed  design  approach  has  proven  to  be  feasible  in  the  experiment  with  the 
book  simulator,  it  still  needs  to  be  tested  in  the  real-time  simulator.  Among  all  the  possible 
practical  issues,  we  foresee  three  that  may  need  to  be  addressed  further.  First,  the  perform¬ 
ance  of  the  designed  controller  may  need  to  be  improved.  As  described  in  Section  5.2,  the 
baseline  controller  is  designed  such  that  the  resulting  stability  region  needs  to  be  maximized. 
While  this  would  provide  the  largest  region  in  the  state  space  for  the  upgraded  controller  to 
explore  new  functionality,  it  may  also  make  the  convergence  to  the  glide  path  slow  when  the 
baseline  control  is  in  charge.  If  this  is  the  case,  constraints  on  the  pole  location  need  to  be 
imposed  to  enforce  the  minimum  decay  rate.  On  the  other  hand,  merely  restricting  the  pole 
location  for  a  minimum  decay  rate  may  cause  the  closed-loop  system  to  oscillate.  When  this 
is  not  acceptable,  the  magnitudes  of  the  imaginary  part  of  the  poles  must  be  limited  as  well. 
All  of  these  can  be  tested  with  different  specifications  of  pole  location. 

Second,  the  accuracy  of  the  identified  model  may  need  to  be  further  improved.  Since  the 
identified  model  is  used  to  represent  the  actual  system,  it  is  always  understood  that  some  dy¬ 
namics  of  the  actual  system  may  not  be  included  in  the  identified  model.  Therefore,  the  actual 
system  can  be  considered  as  the  combination  of  the  identified  model  and  some  uncertain  dy¬ 
namics.  While  the  state  feedback  control  is  usually  robust  in  the  sense  of  controlling  the  sys¬ 
tem  with  uncertainties,  the  stability  region  built  upon  the  identified  model  may  be  problem¬ 
atic.  Furthermore,  since  the  internal  states,  Se, ,  SaiI ,  and  5rdr ,  may  not  be  measurable,  an 

observer  may  be  needed  to  estimate  these  states.  The  accuracy  of  the  observed  states,  of 
course,  will  depend  on  the  precision  of  the  model.  Therefore,  the  accuracy  of  the  identified 
model  will  need  to  be  improved  if  it  gives  poor  state  estimates  or  false  results  on  the  safety 
check  against  the  stability  region.  (See  [Seto  99]  for  information  about  the  safety  check.) 
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Finally,  the  non-linearities  in  the  real-time  simulator  may  cause  problems  in  the  safety  check 
and  system  performance  evaluation.  Although  the  non-linearities  will  result  in  model  inaccu¬ 
racy  as  mentioned  above,  they  may  not  be  addressed  by  improving  the  identified  model.  In 
this  case,  one  needs  to  consider  shrinking  the  stability  region  to  make  it  more  conservative 
and  to  experiment  with  the  performance.  If  the  model  variation  caused  by  the  non-linearities 
can  be  characterized  by  a  bounded  term,  robust  control  methodology  may  then  be  applied. 

In  summary,  the  designed  controller  and  the  derived  stability  region  should  work  well  in  the 
real-time  simulator  after  some  second-order  effect  is  addressed. 
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Appendix  A:  Relationship  Among 

Coordinate  Frames 


Coordinate  Transformation  (*,  y)  —>  (X,Y) : 

Given  a  point  ( xp,yp )  in  the  old  coordinate  frame  (x,y) ,  find  its  expression  (X  P,YP)  in  a 
new  coordinate  frame  (X,7). 

X p  =  xp  cos6  +  yp  sin# 

Yp  =-xp  sin#  + cos# 

or 

Xp  cos#  sin#  xp 

Yp  -sin#  cos#  y  p 

Mapping  from  ECI  to  NED  RNE : 

Celestial  longitude  angle  l :  the  angle  measured  in  the  ECI  y-z  plane  (Earth’s  equatorial 
plane),  from  the  negative  z-axis  in  a  counterclockwise  direction  when  viewed  from  the  North 
Pole. 

Geodetic  latitude  angle,  // :  the  angle  between  the  line  that  is  normal  to  the  Earth’s  surface 
and  passes  the  point  of  interest,  and  the  ECI  y-z  plane. 

These  two  angles  are  depicted  in  Figure  31 . 
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Figure  31:  Celestial  Longitude  Angle  and  Geodetic  Latitude  Angle 
Let  qNED  and  qECI  be  a  vector  expressed  in  the  NED  frame  and  in  the  ECI  frame,  respec¬ 
tively.  Then,qNED  =  RNEq  ECI , 

Rne  =  Rfj  Ri 
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Mapping  from  NED  to  ABC  RAN  : 

Starting  from  NED, 

1 .  Rotate  \ff  about  the  z-axis:  yawing  with  nose  right  positive. 

2.  Rotate  6  about  the  new  y-axis:  pitching  with  nose  up  positive. 

3.  Rotate  (p  about  the  new  x-axis:  rolling  with  right  wing  down  positive,  where  <p,  6  and  if/ 
are  called  the  Euler  angles.  Let  p  ABC  and  pNED  be  a  vector  expressed  in  the  ABC  and 
NED  frames,  respectively.  Then,  pABC  =  RawPnf.d  > 


‘i 

0 

0 

cos  # 

0 

-sin# 
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sin  y 
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cos  <p 

sin  (j) 
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-  cos  0  sin  yr  +  sin  0  sin  #  cos  yr  cos  0  cos  yr  +  sin  0  sin  #  sin  yr  sin  <f>  cos  6 
sin  0  sin  y/ +  cos  0  sin#  cosy/-  -  sin  <f>  cos  \ff  +  cos  <p  sin  8  sin  y/  cos^cos# 
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Mapping  from  ABC  fixed-body  axes  to  ABC-Wind  axes  Rm : 


Starting  from  the  ABC  fixed-body  frame, 

1.  Rotate  about  the  y-axis  for  the  angle  of  a  . 

2.  Rotate  about  the  new  z-axis  for  the  angle  of  p  . 

Let  pWIND  and  PBOdy  be  a  vector  expressed  in  the  ABC-Wind  and  ABC  fixed-body  frames, 
respectively.  Then  pWIND  =  RWBpBom , 

R\VB  =  Rp  R a 


COS  p 

sin  P 

o' 

cos  a 

0 

sin  or 

cos  or  cos  P  sin  P  sin  or  cos  P 

-sin  P 

cos  p 

0 

0 

1 

0 

= 

-cos  or  sin/?  cos  P  -sin  or  sin 

0 

0 

1 

-  sin  or 
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cos  or 

-sin  or  0  cos  or 
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Appendix  B:  Derivation  of  Equations  of 

Motion  of  an  Aircraft 


In  this  section,  we  derive  the  equations  of  motion  of  an  aircraft.  Since  an  aircraft  can  be 
treated  as  a  rigid  body  moving  in  three-dimensional  space,  we  first  review  some  materials  on 
the  motion  of  a  rigid  body,  and  then  derive  the  equations  of  motion  for  an  aircraft  in  the  con¬ 
text  of  a  rigid  body. 

B, 1  Motion  of  a  Rigid  Body 

Motion  between  coordinate  systems:  Let  A  and  B  be  oriented  Euclidean  spaces, 

q,  c  €  A,  and  q  B  e  B  .  A  motion  of  B  relative  to  A  is  a  mapping  smoothly  depending  on  t: 

Dt  ;  B  — )  A, 

which  preserves  the  metric  and  the  orientation.  Every  motion  D,  can  be  uniquely  written  as 
the  composition  of  a  rotation  Rt  :B  A  and  a  translation  C,  :  A  — »  A ,  where 

C, q(f)  =  q(f)  +c(f) .  Namely, 

D, =C,R,  =>  q(0  =  Aqfi(0  =  tf,qfi(t)  +  c(0- 

A  motion  D,  is  called  a  rotation  if  it  takes  the  origin  of  B  to  the  origin  of  A{  i.e.,  if  D,  is  a  lin¬ 
ear  operator).  A  motion  D,is  called  translational  if  the  mapping  does  not  depend  on  t  (i.e., 

R,  =  R0=  R  andZ),qB  =  RqB  +c(t) ). 

Motion  of  a  vector  in  different  coordinate  systems: 

q(t)  =  R,qB(t)+e(f)  =>  q(t)  =  R,qB(t)  +  R,qB(t)+c(t) 

RtqB  =  cox(q-c)  =  <DX(  R,qB) :  the  transferred  velocity  of  rotation 
R,qB  (t) :  the  relative  velocity 

c(t) :  the  velocity  of  the  origin  of  the  moving  coordinate  system 

Where  to  is  the  angular  velocity  of  frame  B  with  respect  to  frame  A  expressed  in  frame  A, 
and  is  derived  from  a  skew-symmetric  matrix  Q.  as 
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0 


-co. 

0 

cor 


Q=R,Rj  = 


CO 

z 

—  CO.. 


COy 

-*>x 

0 


£2q  =  toxq,  Vqe  R 3 


Suppose  c (t)  =  0,  Vr  >  0  .  Let  coB  be  the  angular  velocity  of  frame  B  with  respect  to  frame  A 
expressed  in  frame  B.  Then  ro  =  R,(o  B,  and  the  transferred  velocity  can  be  further  written  as 

R,qB  =  G>xq  =  (fl,coB)x(/?,qB)  =  /?,(0BxqB. 

Define  a  skew-symmetric  matrix  QB  such  that  WjXv  =  fiBv,  Vv  e  R} .  Then 


RiQb  ~  ^  ^ b  ~  Rf 

B.2  Derivation  of  Aircraft  Equations  of  Motion 

Expression  of  a  velocity  vector  in  a  different  coordinate  frame:  Let  A  and  B  be  two  coor¬ 
dinate  frames  and  R  be  the  rotation  matrix  mapping  A—>B.  Then  we  may  define  a  velocity 
vector  as  follows: 

vA  with  respect  to  frame  A  expressed  in  A 
\ba  with  respect  to  frame  A  expressed  in  B 
vB  with  respect  to  frame  B  expressed  in  B 

V*=P>1.  v®=/?vA,  VB  =^-(RpA)  =  (oxpB+RpA 

at 

where  to  is  the  angular  velocity  of  frame  A  with  respect  to  frame  B,  and  pA  and  pB  are  the 
expressions  of  the  position  vector  in  frame  A  and  frame  B,  respectively. 

Example:  Let  frame  A  and  frame  B  be  aligned  along  the  z  axes  as  shown  in  Figure  32  (posi¬ 
tive  z  points  outwards).  Suppose  A  is  stationary  and  B  rotates  about  the  z-axis  at  an  angular 
velocity  6 .  Consider  a  vector  fixed  on  the  X-axis  of  frame  B  with  length  r  and  expressions 
pA  and  pB  in  A  and  B,  respectively.  Then, 
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sin# 
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-sin# 

cos# 
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0 

1 

0 

0_ 

where  R:A-*B.  Hence  the  velocities  can  be  derived  as  the  following: 
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Figure  32:  Vector  Conversion  Between  Frames 

External  forces  in  the  equations 

Propulsion 
Aerodynamic  forces 
Gravitational  attraction  (gravity) 

State  variables  in  the  equations 

Position  vector  in  the  inertial  frame  (ECI):  p  =  [ px ,  py ,  pz]T  . 

Relative  velocity  of  aircraft  eg  with  respect  to  air  mass  and  expressed  in  ABC  frame: 

VB  =K,  VyB,VZBf  . 

Absolute  angular  velocity  of  ABC  frame  expressed  in  ABC:  coB  =  [CQXB ,  0)'B,  6)B]T . 
Attitude  variables  -  Euler  angles  <J>  =  [(p,  6,  lff]T  . 

Dynamic  Equation-Translational  Motion 
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Let  F  be  the  combined  propulsion  and  aerodynamic  forces  and  p  be  the  position  vector  of  the 
aircraft  center  of  gravity,  both  expressed  in  the  ECI  coordinate  frame.  Then  Newton’s  second 
law  implies  that 

F  +  mg  =  — (mp)  =>  F  +  mg  =  mp  +  mp 
dt 

For  aircraft,  we  ignore  the  change  of  mass  (i.e.,  m  =  0 ).  To  derive  an  expression  for  the  vec¬ 
tor  p  ,  let  RECI  be  yet  another  coordinate  system,  which  is  centered  at  the  origin  of  the  ECI 

and  coincident  with  ECI  at  time  t  =  0,  and  rotates  with  the  Earth  about  its  x  axis.  Further¬ 
more,  suppose  Rer  and  Rra  are  the  rotation  matrices  that  result  in  the  mapping 
Rer  :RECI  ->  ECI  and  Rra  :  ABC  -4  RECI .  Let  R  be  the  rotation  matrix  mapping: 

ECI  -»  ABC  (i.e.,  R  =  RANRNE ).  It  follows  that  RT  =  R^R^  ,  and  we  have 

P  =  ^erPr  ^  P  =  ^erPr  ^erPr 

where  pB  is  the  position  vector  expressed  in  RECI  frame.  Clearly,  the  first  term  Rer\>r  is 

the  transferred  velocity  of  the  RECI’s  rotation  (see  Section  B.l  for  the  definition  of  trans¬ 
ferred  velocity),  and  it  can  be  written  as  Rer p  R  =  to  E  x  p  with  to  E  the  absolute  angular  ve¬ 
locity  of  the  Earth’s  rotation  expressed  in  the  ECI  frame.  The  other  term  Rer pB  is  the  rela¬ 
tive  velocity  that  can  be  expressed  as  Rer pr  =  Rc/;(Rr.,vb)  =  RTvB  ■  Finally  the  equation  of 
p  can  be  written  as  follows: 

p  =  to£xp  +  RrvB 

Differentiating  p ,  we  obtain  the  following: 
p  =  a>£xp  +  to£xp  +  /?rvB  +Rt\b 

Since  the  Earth  is  rotating  at  a  constant  angular  velocity,  <b£  =  0 .  We  also  realize  that  the 
third  term,  RT  v  B ,  is  the  transferred  acceleration  due  to  the  rotation  of  the  ABC  frame. 
Hence,  RT\B  =(RT(aB)x(RT\B)  =  RT (i oB  xvB) ,  where  RT coB  and  RT\B  are  the  ABC 

angular  velocity  and  aircraft  cg’s  relative  velocity  with  respect  to  the  air,  both  expressed  in 
the  ECI  coordinate  frame.  Then  the  equation  of  p  becomes 

p  =  co£x(coBxp  +  /?7'vB)  +  Rr(toBxvB)  +  RrvB 

and  the  equation  of  Newton’s  second  law  can  be  written  as 

F  +  mg  =  m[©£  x (coB  xp  +  Rr vB )  +  RT (toB  x  vB)  +  RT vB ] 
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Let  FB  be  the  force  F  expressed  in  ABC  coordinate  frame(i.e.,  FB  =  RF ).  Then  the  above 
dynamic  equation  can  be  written  in  ABC  by  pre-multiplying  matrix  R  as 

FB+Rmg  =  m[R(aEx(<aBxp  +  RT\B)  +  (oBx\B+yB] 

=  m[R(»Ex((oB  xp)  +  (R(oE)x\ B  +  toB  xvB  +  vB] 

Finally,  we  obtain 

P 

vB  =-J--(coB+/?co£)xvfi+/?[g-co£x(coBxp)] 
m 


Dynamic  Equation-Angular  Motion 

Let  M  be  the  angular  momentum  and  T  be  the  total  torque  acting  about  the  aircraft  eg,  both 
expressed  in  the  ECI  coordinate  frame.  Then  Newton’s  second  law  implies 


M  =  T 


Since  the  angular  momentum  of  the  aircraft  in  the  ABC  coordinate  can  be  expressed  as  Jo  B , 
where  J  is  the  moment  of  inertial  of  the  aircraft,  then  M  =  RT  Ji o  B  and 
M  =  RT  J(aB  +  RtJ6)b  .  Again,  the  first  term  is  recognized  as  the  transferred  momentum  and 
it  can  be  written  as  RTJ(oB  =  (RTtnB)xM,  =  RT ( toB  x(Ja>B)) .  Then  the  dynamic  equation 
can  be  rewritten  as  follows: 

T  =  RT (i db  x(7<ob))+  RtJ6)b 

Let  Tb  be  the  ABC  expression  of  T  (i.e.,  TB  =  RT ).  Then  the  dynamic  equation  can  be  ex¬ 
pressed  in  ABC  coordinate  as  follows: 

d>B  =-/-‘[<oBx(ya)B)]+y-ITB 

Kinematics  Equations 


Three-variable  attitude  equation: 

Define  the  Euler  angles  [(j),6,y/]  between  the  ECI  frame  and  the  ABC  frame,  and  let  RAE  be 
the  rotation  matrix  mapping  from  ECI  to  ABC.  Then, 
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and 


&B  ~  ^AE^AE  ~~^AE^AE  ^  ^AE~  &B^AE 

where  QB  is  the  skew-symmetric  matrix  such  that  toB  xq  =  £2Bq  .  By  evaluating  the  ele¬ 
ments  (1,2),  (1,3),  and  (2,3)  in  the  above  equation,  we  obtain  the  following: 
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tan  0  sin  <p 

tan  #cos  (j> 

coB 

0 

= 

0 

COS  Cf) 

-sin  0 
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Four-variable  attitude  equation: 


The  Euler  angles  can  be  replaced  by  the  so-called  quaternion  four-variable  representation, 
<1  =  [<7o  >  9i  >  92  >  93  ]T  •  Then  the  rotation  matrix  RAE  from  the  ECI  frame  to  the  ABC  frame 
can  be  expressed  in  terms  of  q  as  follows: 


Rae  ~ 


2  2  2  2 
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and  the  relation  between  q  and  <D  is  given  by  the  following: 

0  =  tan-1— ,  0  =  -sin-'  fc13,  ^  =  tan_1-^l 

^33  1 

with  by  is  the  element  in  the  ith  row  and  )th  column  of  RAE  ,  and 

q0  =  ±[cos(0  /  2)  cos(0  /  2)  cos(^>/  2)  +  sin(0  /  2)  sin(0  /  2)  sin($?  /  2)] 
g,  =  ±[sin(^  /  2)  cos(0  /  2)  cos(^o  /  2)  -  cos(0  /  2)  sin(0  /  2)  sin(#>/  2)] 
g2  =  ±[cos(0  /  2)  sin(0  /  2)  cos(q?  /  2)  +  sin(0  /  2)  cos(0  /  2)  sin(^>  /  2)] 

<j3  =  ±[cos(0  /  2)  cos(0  /  2)  sin(^  /  2)  -  sin(0  /  2)  sin(0  /  2)  cos(^)  /  2)] 

Therefore,  the  kinematics  equations  in  terms  of  q  are  given  by  the  following: 
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B.3  Summary  of  Modelling 

The  Round-Earth  Equations  in  ECI  are  as  follows: 


P  -  ^£P  +  ^A£V£ 

\B  =  “  RaE^eP  ”  (^B  +  ^AE^E^AE^B  ^A£§(P)  + 

co5  =  —  /  +7  ]Tb 


where 
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Ci\  is  the  x-component  of  the  absolute  angular  velocity  of  the  Earth  in  ECI,  FB  and  T/f  are  the 
resultant  force  and  torque  of  propulsion  and  aerodynamic  forces,  expressed  in  the  ABC 
frame,  and  g(p)  is  the  gravitational  attraction. 


The  Flat-Earth  Equations  in  NED 

In  this  case,  we  also  ignore  the  Earth’s  rate  (i.e.,  to E  =  0 ).  Then  the  equations  become 


Pned  ~  Ranvb 

F 

\B  =  —  QbVB  +  ^AvSo  "* 

m 

to B  —  —J  ^SiBJoB  +J  *T b 

6  =  T(0)tofi 


where 


1  tan  0  sin  (j)  tan  0  cos  0 

0 

~Pn] 

¥(<!>)  = 

0  cos  (j)  -sin$> 

0  sin^/cos#  cos0/cos0 

?  g0 

0 

32.17  ft/s2_ 

*  Pned  ~ 

- 1 

_ 1 

and  the  Euler  angles  are  defined  between  the  NED  frame  and  the  ABC  frame. 
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Appendix  C:  Data  Collection  /  Protection 

Methodology 


The  legacy  F-16  manned  virtual  desktop  simulation  mentioned  in  the  Introduction  is 
an  unclassified,  but  International  Tariffs  in  Arms  Regulations  (ITARS)-restricted 
item.  This  meant  that  non-U.S.  nationals  could  not  have  access  to  the  artifact,  which 
presented  some  difficulties  in  the  collection  and  analysis  of  the  data.  This  appendix 
outlines  the  procedures  that  the  team  employed  to  protect  the  information  represented 
in  the  simulation  artifact. 

The  simulation  is  a  reduced  order  F-16  aircraft  model  running  in  a  soft  real-time 
desktop  computational  environment.  This  simulation  is  assumed  by  the  authors  to  be 
representative  of  F-16  flight  characteristics  (indeed,  of  the  entire  class  of  fighter  air¬ 
craft),  but  it  should  not  be  considered  a  ground  truth  information  source  for  the  F-16 
aircraft.  The  simulation  was  not  certified  as  correct  by  any  process,  and  no  flight  test 
evaluation  by  any  qualified  pilot  was  performed.  Any  relationship  between  data  col¬ 
lected  from  this  simulation  and  real-world  ground  truth  is  therefore  unknown. 

In  addition,  the  flight  envelope  inspected  by  this  process  has  no  tactical  significance 
to  F-16  flight  operations.  All  data  were  collected  in  a  “gear  down”  landing  configura¬ 
tion.  The  simulation  was  initialized  with  the  aircraft  within  3  degrees  laterally  of  the 
runway  localizer,  at  1120  feet  AGL  altitude  (3400  feet  MSL  altitude),  and  a  velocity 
of  218  KCAS.2  The  aircraft  was  hand-flown  by  engineering  personnel  from  this  point 
to  touchdown. 

Data  were  extracted  from  the  model  at  a  50  Hz  rate.  The  extracted  data  included  the 
horizontal  and  vertical  glide-path  deviations,  aircraft  Euler  angles  and  rates,  altitude, 
airspeed,  vertical  velocity,  angle-of-attack,  pilot  control  inputs  and  control  surface 
positions. 

A  U.S.  national  was  the  identifier,  operator,  and  data  collector  for  this  system.  At  no 
time  did  non-U.S.  nationals  have  unrestricted  access  to  either  the  source  code  or  ex¬ 
ecutable  simulation  artifacts.  Once  collected,  the  data  were  transferred  to  other  com¬ 
puters  for  analysis  and  model  identification. 


2  The  acronyms  AGL,  MSL,  and  KCAS  are  defined  as  follows:  AGL  is  above  ground  level;  MSL  is 
mean  sea  level;  and  KCAS  is  knots  calibrated  air  speed. 
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